##================================================## ## ## ## Computational Social Science ## ## Doctoral Seminar (16:194:630) ## ## Spring 2026, Rutgers University ## ## Katya Ognyanova, katya@ognyanova.net ## ## ## ##================================================## # ================ Machine Learning in R ================ # Packages we would need today: # (remove # and run this if they are not already installed) # install.packages("tidymodels") # Machine learning workflow # install.packages("randomForest") # Random Forest classifier # install.packages("kknn") # KNN classifier library(tidymodels) library(ggplot2) # ================ ~~ Examine the data ================ # Today we'll use a dataset describing penguins who live on the Palmer Archipelago: data(penguins) penguins # Rename variables to make life easier: colnames(penguins) <- c("species", "island", "bill_length", "bill_depth", "flipper", "mass", "sex") # Look at the data! View(penguins) # Look at the penguins! browseURL("https://www.bas.ac.uk/wp-content/uploads/2015/04/Penguin-heights-736x419.jpg") browseURL("https://allisonhorst.github.io/palmerpenguins/reference/figures/culmen_depth.png") # Check the size of our data: dim(penguins) # Check how many missing values we have in each variable: colSums(is.na(penguins)) # Quick summary of all variables: summary(penguins) # Check correlations among numeric variables: penguin_cor <- cor(penguins[, 3:6], use = "complete.obs") penguin_cor # Plot mass by flipper length, colored by species ggplot(penguins, aes(x=flipper, y=mass, color = species)) + geom_point(size=2, alpha=0.7) + scale_color_manual(values = c("darkorange","darkred","skyblue3")) + theme_bw() # Plot mass by flipper length, split by species, colored by sex ggplot(penguins, aes(x = flipper, y = mass)) + geom_point(aes(color = sex),size=2, alpha=0.7) + scale_color_manual(values = c("darkred", "deepskyblue4")) + facet_wrap(~species)+ theme_bw() # Plot bill length by bill depth, colored by species ggplot(penguins, aes(x=bill_length, y=bill_depth, color = species)) + geom_point(size=2, alpha=0.7) + scale_color_manual(values = c("darkorange","darkred","skyblue3")) + theme_bw() # Let's create a "long" version of our dataset for plotting: # We are combining columns 3 to 6, creating two new columns # The old variable names are now in "measure", the values in "score" penguin_long <- pivot_longer(penguins, cols = 3:6, names_to = "measure", values_to = "score") penguin_long ggplot(penguin_long, aes(x=score, fill = species)) + geom_histogram(alpha = 0.7, bins=20, position = "identity") + scale_fill_manual(values = c("darkorange","darkred","skyblue4")) + facet_wrap(~ measure, scales = "free_x") + theme_bw() # ================ ~~ Our penguin ML workflow ================ # tidymodels gives us nice consistent steps to follow for machine learning: # # 1. SPLIT: Divide data into training and test sets before doing anything else. # The test set is left aside until final evaluation. # initial_split() --> training() / testing() # # 2. RECIPE: Define the outcome, the predictors, and any preprocessing steps # (imputation, scaling, encoding, etc.). The recipe is learned # from training data only -- this prevents information from the # that test set we left aside from leaking into our model. # recipe() --> step_() --> prep() --> bake() # # 3. MODEL: Specify what kind of model to use, what it should do # (classification or regression), and which package runs it. # rand_forest() / nearest_neighbor() / etc. # --> set_mode() --> set_engine() # # 4. WORKFLOW: Bundle the recipe and model together so preprocessing and # fitting are always applied consistently as a single unit. # workflow() --> add_recipe() --> add_model() --> fit() # # 5. EVALUATE: Generate predictions on the test set and assess performance # using metrics appropriate to the task. # predict() --> conf_mat() --> summary() # # 6. TUNE: Use cross-validation to try different hyperparameter values # and select the best one -- all within the training data. # vfold_cv() --> tune_grid() --> select_best() --> finalize_workflow() # --> last_fit() to train final model and evaluate on test set once. # # The key thing to keep in mind all the time: test data is used only once, at the end. # All modeling decisions (preprocessing, model choice, hyperparameters) will be made # using the training data only. # ================ ~~ Preprocess the data ================ # Because splitting the data is random, we set a seed to make sure # you do the same split and get the exact same results as I do set.seed(111) # Below, we are doing a simple random splitting of our dataset, with # 70% of the data used for training and 30% withheld to test the model. # If we had a large class imbalance (lots of cases in some categories # and very few in others), we could do stratified sampling, which would # do the training/testing data split within each category and combine. # For example, we could include "strata=species" below to split within species. penguins_split <- initial_split(penguins, prop = 0.7) penguins_split # Get our training and testing data based on the split: penguins_train_data <- training(penguins_split) penguins_test_data <- testing(penguins_split) # Preprocess the data for classification. # The tidymodels package offers a way to do that using a 'recipe'. # It specifies our outcome variable, predictor variables, as well as # what we want to do with the data to filter and clean it for analysis. # We can then "bake" different datasets using this recipe. # First, we will specify the outcome variable and the predictors (features) we want to use. # We are using a formula, much like we did when we ran regression models. # Specifying what the outcome and predictors are going to be is important: # the recipe needs to know in order to apply cleaning steps correctly. # For example, we only want to apply some pre-processing steps (e.g. dropping # collinear variables) to the predictors, and never to the outcome. penguins_recipe <- recipe(species ~ bill_length + bill_depth + flipper + mass + sex, data = penguins_train_data) # Next, we'll add some data cleaning steps to our initial recipe. # Replace missing values in the data with imputed values from other observations: penguins_recipe <- step_impute_median(penguins_recipe, all_numeric_predictors()) penguins_recipe <- step_impute_mode(penguins_recipe, all_nominal_predictors()) # Remove any variables that are too highly correlated. Apply to all numeric predictors. penguins_recipe <- step_corr(penguins_recipe, all_numeric_predictors()) # Normalize the data so it has a mean of 0. Apply to numeric predictors. penguins_recipe <- step_center(penguins_recipe, all_numeric_predictors()) # Scale the data to have a standard deviation of 1. Apply to numeric predictors. penguins_recipe <- step_scale(penguins_recipe, all_numeric_predictors()) # Note that there are many other recipe steps we could add based on our data and goals. # Some useful examples include: # # step_naomit: drop rows with missing values instead of imputing # step_impute_knn: impute using k-nearest neighbors (vs. median/mode) # step_nzv: remove near-zero variance predictors # step_pca: reduce dimensionality via principal components # step_log/step_sqrt: transform skewed numeric variables # step_poly: add polynomial terms for numeric predictors # step_interact: add interaction terms between predictors # step_upsample/step_downsample: address class imbalance in the outcome # # Steps can also be applied to specific subsets of variables using selectors: # all_predictors(), all_outcomes(), all_numeric_predictors(), all_nominal_predictors() # or by name using starts_with(), contains(), matches(), etc. # Calculate everything we need to execute the recipe (e.g. variable means, etc.) penguins_recipe_prepped <- prep(penguins_recipe) # Our recipe is all done! penguins_recipe_prepped # Apply the prepped recipe to the training data. # You can use a shorter equivalent function: juice(penguins_recipe_prepped) penguins_training <- bake(penguins_recipe_prepped, new_data=NULL) # Apply the same prepped recipe to our test data: penguins_testing <- bake(penguins_recipe_prepped, new_data=penguins_test_data) # The new datasets are just our original data with the additional # steps described in the recipe applied to them. penguins_training penguins_testing # ================ ~~ Training model ================ # There are many classification techniques (e.g. linear and logistic regression, # classification trees, random forests, k-nearest neighbors, support vector machines, # neural networks, and others). We will not cover any of them in much detail today. # Here, we will try two classifiers: random forest (RF) and k-nearest neighbors (KNN). # Specify the RF model: penguins_model_rf <- rand_forest(trees = 100) # What model? penguins_model_rf <- set_mode(penguins_model_rf, "classification") # What should it do? penguins_model_rf <- set_engine(penguins_model_rf, "randomForest") # Which package should run it? # Train the model on the data: penguins_fit_rf <- fit(penguins_model_rf, species ~ ., data = penguins_training) penguins_fit_rf # Remember, in the formula above, species ~ . tells R to use as a predictor # every available variable in the data other than the outcome (species). # This saves us some time typing the full formula. # Specify the KNN model: penguins_model_knn <- nearest_neighbor(neighbors = 5) penguins_model_knn <- set_mode(penguins_model_knn, "classification") penguins_model_knn <- set_engine(penguins_model_knn, "kknn") # Train the model on the data: penguins_fit_knn <- fit(penguins_model_knn, species ~ ., data = penguins_training) penguins_fit_knn # Notice that one thing we do not do here is examine how the model results # change for different parameter values like trees or neighbors above. # Typically we would want to tune the parameters to get the best # possible model we can generate. We'll show that a bit further down. # ================ ~~ Data prediction ================ # Predict the penguin species using the trained model and training data: penguins_predict_train_rf <- predict(penguins_fit_rf, penguins_training) penguins_predict_train_knn <- predict(penguins_fit_knn, penguins_training) # Predict the penguin species using the trained model and test data: penguins_predict_test_rf <- predict(penguins_fit_rf, penguins_testing) penguins_predict_test_knn <- predict(penguins_fit_knn, penguins_testing) # Add the predictions to our training dataset: penguins_training$predict_rf <- penguins_predict_train_rf[[1]] penguins_training$predict_knn <- penguins_predict_train_knn[[1]] penguins_training # Add the predictions to our test dataset (more interesting!): penguins_testing$predict_rf <- penguins_predict_test_rf[[1]] penguins_testing$predict_knn <- penguins_predict_test_knn[[1]] penguins_testing # We can also obtain predicted probabilities for each possible outcome: penguins_prob_rf <- predict(penguins_fit_rf, penguins_testing, type="prob") penguins_prob_rf penguins_prob_knn <- predict(penguins_fit_knn, penguins_testing, type="prob") penguins_prob_knn # We can add these to our testing dataset as well: names(penguins_prob_rf) <- c("adelie_rf", "chinstrap_rf", "gentoo_rf") penguins_testing <- cbind(penguins_testing, penguins_prob_rf) names(penguins_prob_knn) <- c("adelie_knn", "chinstrap_knn", "gentoo_knn") penguins_testing <- cbind(penguins_testing, penguins_prob_knn) # ================ ~~ Data validation ================ # A "confusion matrix" for our data shows us the actual and predicted values. # Each column includes the cases that actually belong to a category. # Each row shows what category the model predicted for those cases. # (note that in some textbooks, rows and columns are reversed) # Not surprisingly, the classifiers perfectly predict their training data: conf_mat(penguins_training, truth = species, estimate = predict_rf) conf_mat(penguins_training, truth = species, estimate = predict_knn) # What we are interested in is how well they do on the withheld test data: penguin_cm_rf <- conf_mat(penguins_testing, truth = species, estimate = predict_rf) penguin_cm_knn <- conf_mat(penguins_testing, truth = species, estimate = predict_knn) penguin_cm_rf penguin_cm_knn # From the confusion matrix, we can obtain a number of model fit metrics: summary(penguin_cm_rf) summary(penguin_cm_knn) # We can see that KNN has 100% accuracy on test data compared to 97% for RF. # Some of the metrics we might be interested in include: # # accuracy: # Percent of the cases that were predicted correctly. # E.g. what percent of penguins were sorted in the right category? # # precision: # Correctly predicted in category / Total predicted in category (correctly or not). # E.g. what percent of predicted Adelie penguins are really Adelie penguins? # # sensitivity/recall: # Correctly predicted in category / Total truly in that category (predicted or not). # E.g. what percent of real Adelie penguins were predicted to be Adelie penguins? # # f1-score/f-measure: harmonic mean of precision and recall. # The score tries to balance precision and recall. # # specificity: Correctly predicted NOT in category / Total truly NOT in category. # E.g. what percent of penguins we said were not Adelie were really not Adelie? # # Kappa: accuracy over and above the accuracy we can expect by chance. # Values closer to 1 are preferable. # # # Depending on what we are predicting, we may prioritize different metrics. # For example, if we are categorizing sick people to treat, we may be more # concerned about recall and less about precision (more concerned that we # identify everyone who is really sick, and less worried about false positives). # ================ ~~ Workflow ================ # The tidymodels package allows us to create a workflow that # combines our recipe and the training model. # We can then use that workflow with different datasets. # Note that workflows take the unprepped recipe and handle prep internally, # which ensures preprocessing is applied correctly within each fit call. penguins_workflow <- workflow() penguins_workflow <- add_model(penguins_workflow, penguins_model_knn) penguins_workflow <- add_recipe(penguins_workflow, penguins_recipe) # Fit the workflow on the training data (this is what we would normally do): penguins_work_fit <- fit(penguins_workflow, penguins_train_data) penguins_work_fit # Here we fit the same workflow on the test dataset, just to demonstrate # that workflows can be reused. In practice, we would NOT train on test data. penguins_work_fit <- fit(penguins_workflow, penguins_test_data) penguins_work_fit # ================ ~~ Cross-validation ================ # Cross-validation involves splitting our data into training and testing # datasets multiple times, then training and validating the model on each split. # K-fold (or v-fold) cross-validation involves splitting the data into # k equal parts called "folds". We then use k-1 of the folds as our # training data, and the remaining one fold as the test data. # We repeat that k times, independently training and testing the model. # Our validation results are aggregated across all k validations. # We can use k-fold cross-validation to tune model parameters, # and/or to validate the model. # People often use k=5 or k=10 in cross-validation. # Cross-validation with 10-fold split: penguins_val_split <- vfold_cv(penguins, v = 10) penguins_val_split # We can get each of the 10 splits separately if we want to: penguins_val_split$splits[[1]] # Create a workflow bundling the RF model with our recipe. # Using a workflow here ensures the recipe preprocessing steps are applied # consistently within each fold, rather than leaking information across folds. penguins_workflow_rf <- workflow() penguins_workflow_rf <- add_model(penguins_workflow_rf, penguins_model_rf) penguins_workflow_rf <- add_recipe(penguins_workflow_rf, penguins_recipe) # Fit the Random Forest workflow for each of the 10 data splits: # The save_pred parameter makes sure we save the predictions from each fold. penguins_val_fit <- fit_resamples(penguins_workflow_rf, resamples = penguins_val_split, control = control_resamples(save_pred = T)) penguins_val_fit # Get predicted values: penguins_val_predict <- collect_predictions(penguins_val_fit) penguins_val_predict # Get the confusion matrix (averaged across resamples): penguins_val_cm <- conf_mat_resampled(penguins_val_fit, tidy=F) penguins_val_cm # Get metrics from the averaged confusion matrix: summary(penguins_val_cm) # ================ ~~ Parameter tuning ================ # Earlier, we user a KNN model with 5 neighbors. # But, how do we know that neighbors=5 is the right choice for KNN? # We can 'tune' the parameter by trying multiple values and picking the best one. # Specify a KNN model with neighbors marked for tuning. # Basically we're saying we don't know what value to use and want to try a few. penguins_model_knn_tune <- nearest_neighbor(neighbors = tune()) penguins_model_knn_tune <- set_mode(penguins_model_knn_tune, "classification") penguins_model_knn_tune <- set_engine(penguins_model_knn_tune, "kknn") # Create a workflow with the tunable model: penguins_workflow_tune <- workflow() penguins_workflow_tune <- add_model(penguins_workflow_tune, penguins_model_knn_tune) penguins_workflow_tune <- add_recipe(penguins_workflow_tune, penguins_recipe) # Define the values of neighbors we want to try: knn_grid <- tibble(neighbors = c(1, 3, 5, 7, 10, 15, 20)) # Try each value across the 10 cross-validation folds: # tune_grid takes each data split and tests the corresponding k on it penguins_tune_fit <- tune_grid(penguins_workflow_tune, resamples = penguins_val_split, grid = knn_grid, metrics = metric_set(accuracy, precision, recall)) # Look at performance for each parameter value: collect_metrics(penguins_tune_fit) # See the best value of neighbors for different metrics: show_best(penguins_tune_fit, metric = "accuracy") show_best(penguins_tune_fit, metric = "recall") show_best(penguins_tune_fit, metric = "precision") best_neighbors <- select_best(penguins_tune_fit, metric = "accuracy") # Finalize the workflow using the best parameter value: # This step takes the placeholder tuning values and replaces them with our chosen one. penguins_workflow_final <- finalize_workflow(penguins_workflow_tune, best_neighbors) # Fit the final model on the full training data and evaluate on test data: penguins_final_fit <- last_fit(penguins_workflow_final, penguins_split) collect_metrics(penguins_final_fit) # ================ ~~ Regression model ================ # Since we are already familiar with regression from statistics, # regression machine learning models may be easier to understand. # Instead of classifying data, we are predicting numerical values. # We reload the penguins data here and start fresh for the regression example data(penguins) colnames(penguins) <- c("species", "island", "bill_length", "bill_depth", "flipper", "mass", "sex") # Here is a standard regression example: my_regression <- lm(mass ~ species + flipper + sex + bill_length + bill_depth, data=penguins) summary(my_regression) # Now a machine learning model predicting penguin mass: penguins_split_reg <- initial_split(penguins, prop = 0.7) penguins_split_reg # Get our training and testing data: penguins_train_reg <- training(penguins_split_reg) penguins_test_reg <- testing(penguins_split_reg) # Preprocess the data: penguins_recipe_reg <- recipe(mass ~ species + flipper + sex + bill_length + bill_depth, data = penguins_train_reg) penguins_recipe_reg <- step_impute_median(penguins_recipe_reg, all_numeric_predictors()) penguins_recipe_reg <- step_impute_mode(penguins_recipe_reg, all_nominal_predictors()) penguins_recipe_reg <- step_corr(penguins_recipe_reg, all_numeric_predictors()) penguins_recipe_reg <- step_center(penguins_recipe_reg, all_numeric_predictors()) penguins_recipe_reg <- step_scale(penguins_recipe_reg, all_numeric_predictors()) penguins_recipe_reg <- step_dummy(penguins_recipe_reg, all_nominal_predictors()) penguins_recipe_reg_prepped <- prep(penguins_recipe_reg) # And bake the recipe: penguins_train_reg <- bake(penguins_recipe_reg_prepped, new_data=NULL) penguins_test_reg <- bake(penguins_recipe_reg_prepped, new_data=penguins_test_reg) # We can specify a model that uses linear regression: penguins_model_reg_lm <- linear_reg() penguins_model_reg_lm <- set_mode(penguins_model_reg_lm, "regression") penguins_model_reg_lm <- set_engine(penguins_model_reg_lm, "lm") # Train the LM model on the data: penguins_fit_reg_lm <- fit(penguins_model_reg_lm, mass ~ ., data = penguins_train_reg) penguins_fit_reg_lm # And a model that uses random forest: penguins_model_reg_rf <- rand_forest(trees = 100) penguins_model_reg_rf <- set_mode(penguins_model_reg_rf, "regression") penguins_model_reg_rf <- set_engine(penguins_model_reg_rf, "randomForest") # Train the RF model on the data: penguins_fit_reg_rf <- fit(penguins_model_reg_rf, mass ~ ., data = penguins_train_reg) penguins_fit_reg_rf # Predict penguin mass using the trained models: penguins_predict_reg_lm <- predict(penguins_fit_reg_lm, penguins_test_reg) penguins_predict_reg_rf <- predict(penguins_fit_reg_rf, penguins_test_reg) # Add the predictions to our test dataset: penguins_test_reg$predict_lm <- penguins_predict_reg_lm[[1]] penguins_test_reg$predict_rf <- penguins_predict_reg_rf[[1]] penguins_test_reg # Evaluate how well we did in our predictions: metrics(penguins_test_reg, mass, predict_lm) metrics(penguins_test_reg, mass, predict_rf) # Here metrics include: # R Squared: explained variance, higher is better. # Mean Absolute Error: average error between true and predicted values, lower is better. # Root Mean Squared Error: square root of the mean squared error, lower is better. # ================ ~~ BONUS: Penguin Clustering ================ # K-means clustering groups observations by similarity # without specifically being told what the groups should be. # Let's see if clustering can recover the penguin species. data(penguins) colnames(penguins) <- c("species", "island", "bill_length", "bill_depth", "flipper", "mass", "sex") # Keep only the variables we will cluster on: penguins_cluster <- penguins[, c("species", "bill_length", "bill_depth", "flipper", "mass")] # Remove rows with missing values: penguins_cluster <- na.omit(penguins_cluster) # Scale the numeric features so no single variable dominates # due to differences in units or range: penguins_scaled <- scale(penguins_cluster[, 2:5]) # Run K-means with K=3 (because we know there are 3 species) # Setting parameter nstart=25 runs the algorithm from 25 different # random starting points and keeps the best result -- this helps avoid # getting stuck in a suboptimal solution due to an unlucky starting point: penguin_kmeans <- kmeans(penguins_scaled, centers = 3, nstart = 25) penguin_kmeans # Add cluster assignments back to the data: penguins_cluster$cluster <- factor(penguin_kmeans$cluster) # How do the clusters map onto actual species? # Rows are species, columns are clusters: table(penguins_cluster$species, penguins_cluster$cluster) # Visualize: what the algorithm found vs. what the species actually are ggplot(penguins_cluster, aes(x=flipper, y=bill_length, color=cluster, shape=species)) + geom_point(size=2, alpha=0.8) + scale_color_manual(values = c("darkorange","darkred","skyblue4")) + labs(title="K-means clusters vs. actual species") + theme_bw() # ================ THE END ================